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Abstract 

In lattice gauge theory, there exist field transformations that map the theory to the trivial 
one, where the basic field variables are completely decoupled from one another. Such maps 
can be constructed systematically by integrating certain flow equations in field space. The 
construction is worked out in some detail and it is proposed to combine the Wilson flow 
(which generates approximately trivializing maps for the Wilson gauge action) with the 
HMC simulation algorithm in order to improve the efficiency of lattice QCD simulations. 



1. Introduction 

The Nicolai map transforms interacting supersymmetric theories to non-interacting 
ones [1]. Supersymmetry is considered to be essential for the existence of these field 
transformations in view of the fact that their Jacobian is exactly cancelled by the 
fermion partition function. 

In lattice gauge theory, a natural question to ask is whether there are field trans- 
formations that map the theory to its strong-coupling limit. In particular, if there 
are no matter fields, one is looking for substitutions 

U = F{V) (1.1) 

of the gauge field U in the functional integral whose Jacobian cancels the gauge-field 
action. Similarly to the Nicolai map, this kind of transformation maps the theory to 
a solvable one, but supersymmetry is not required and the Jacobian plays a different 
role. 
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Fig. 1. The proposed simulation algorithm for lattice QCD updates the gauge field 
U in three steps, following the arrows in this diagram, where the Hamilton function 
used in the HMC step has the standard kinetic term and includes the Jacobian of the 
field transformation T (cf. subsect. 2.4). 

On a finite lattice, and if the gauge group is compact and connected, the existence 
of such trivializing maps is guaranteed by a general theorem on volume forms on 
compact manifolds (see ref. [2], Theorem 1.26, for example). One may be inclined to 
assume that these transformations are too complicated to be of any use. However, as 
explained later in this paper, it is possible to build up trivializing maps by integrating 
flows in field space, whose generators satisfy certain partial differential equations. 
The latter are quite tractable and can, to some extent, be solved analytically in the 
pure gauge theory. An application of trivializing maps, which can then be envisaged, 
is the acceleration of lattice QCD simulations. 

The fact that the efficiency of the available simulation algorithms is unpredictable 
has always been a weakness of numerical lattice QCD. Already a while ago, empirical 
studies of the SU(3) gauge theory by Del Debbio, Panagopoulos and Vicari [3] showed 
that the autocorrelation times of observables related to the topological charge of the 
gauge field tend to be large and appear to grow exponentially with the inverse of the 
lattice spacing. Moreover, Schaefer, Sommer and Virotta [4] recently found that the 
situation is, in this respect, essentially unchanged when the sea quarks are included 
in the simulations. 

The rapid slowdown of the simulations at small lattice spacings may conceivably 
be overcome by combining approximately trivializing maps with the HMC simulation 
algorithm [5] (see fig. 1). Since the transformation moves the theory closer to the 
strong-coupling limit, where the HMC algorithm is known to be highly efficient, the 
autocorrelation times are, in general, expected to be reduced in this way. Evidently, 
for the combined algorithm to work out in practice, approximately trivializing maps 
must be found which are fairly simple and programmable. One of the goals of the 
present paper is thus to provide a solution to this problem. 



2. Field transformations 



Most concepts developed in this paper are expected to be widely applicable, but as 
explained above, the case of immediate interest is lattice QCD. In the following, the 
gauge group is therefore taken to be SU(3). Since the quarks will play a spectator 
role, they will be added to the theory only in sect. 6, where the proposed simulation 
algorithm for lattice QCD is discussed. 

2.1 Field space 

The lattice theory is set up on a finite hypercubic lattice A with periodic boundary 
conditions. For notational convenience, the lattice spacing is set to unity. As usual, 
the gauge field variables U (x, fi) € SU(3) are assumed to reside on the links (x, fi) of 
the lattice (where x € A and jj, = 0, . . . ,3). The expectation value of any observable 
0(U) is then given by the functional integral 

(O) = 1 j B[U}0(U)e- s ^ (2.1) 

over the space of gauge fields. In this expression, S(U) denotes the gauge-field action, 
Z the partition function and D[f7] the product of the normalized SU(3)-invariant 
integration measures of the link variables U(x,fi). 

From a purely mathematical point of view, the space of lattice gauge fields is a 
power of SU(3) and therefore a compact connected manifold. Field transformations 
are invertible maps of this manifold onto itself. Such transformations will always be 
required to be differentiable in both directions and orientation-preserving (here and 
below, "differentiable" means "infinitely often continuously differentiable"). 

2.2 Right-invariant differential operators 

Since the link variables U(x,[i) take values in a Lie group, it is natural to express 
differentiations with respect to them through a basis <9" M , a = 1, . . . , 8, of differential 
operators that are invariant under the right-action of the group. The action of these 
operators on a differentiable function f{U) of the gauge field is given by 



{ e tTa U(x,u) if (y, u) = (x,u), 
, U t (3,,u) = \ (2.2) 
t=o [ U(y,u) otherwise, 



where T a are the SU(3) generators (see appendix A). In particular, e?£ transforms 
according to the adjoint representation under the left-action of SU(3). 
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The operators go along with a basis of 1-forms on the field manifold, 

0»„ = -2tr{d[/(x^)C/(x,/.)- 1 T a }, (2.3) 
such that 

d/(C/) = E«V( C/ ) ( 2 - 4 ) 

for all functions /(£/)• 
2.5 Jacobian matrix 

For any given gauge field V(y, v), the transformation (1.1) produces another field 
U (x, n) = [T(y)\(x, fi). When considering such transformations, one needs to distin- 
guish differentiations with respect to V from those with respect U. The associated 
1-forms must also be distinguished. In the following, all symbols carrying a hat 
represent quantities and operations referring to V. 
The Jacobian matrix 

[F*(V)](x,»;y,vr b = -2tr{d b yil/ U(x, v)U(x, /i)" 1 ^} (2.5) 

can be considered to be the kernel of a linear operator acting on link fields with 
values in su(3). In particular, 

e a X! , = Y,^*( v )^,^y,v) ab dl„- (2-6) 

y,v 

Since the functional integration measure D[£7] is proportional to the maximal prod- 
uct of these 1-forms, it follows that 

D[L7] =D[V]detF*(V). (2.7) 

The Jacobian of the map (1.1) is thus det.F*(V). 
If the transformation satisfies 

S{T(V)) - lndet^(y) = constant, (2.8) 

the substitution U — ► V of the integration variables in the functional integral maps 
the theory to the trivial one where the link variables are completely decoupled from 
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one another. The expectation values (2.1) are then given by 



(0) = J D[V](D(T(V)) 



(2.9) 



Such trivializing maps thus contain the entire dynamics of the theory. 

Although the remark is likely to remain an academic one, an intriguing observation 
is that the integral (2.9) can be simulated simply by generating uniformly distributed 
random gauge fields. Subsequent field configurations are uncorrelated in this case 
and there are therefore no autocorrelations in the data series for the observables O. 

2.4 Transformation behaviour of the HMC algorithm 

The HMC algorithm [5] operates on the phase space associated to the field manifold. 
In particular, the transition V — > V in fig. 1 requires the equations of motion derived 
from the Hamilton function 



to be integrated, where tt{x,h) <G su(3) is the canonical momentum of V{x,n). 

Although the complete update algorithm for the field U described by fig. 1 looks 
different, it is in fact equivalent to the HMC algorithm with a non-standard Hamilton 
function. The equivalence can be established by noting that the transformation from 
V to U preserves the symplectic 2-form 



The evolution of the transformed fields is then governed by the Hamilton function 



H{tt,V) = i(7r,7r) + S(F(V)) -lndet^(F) 



(2.10) 




(2.11) 



i.e. f2 = f2, if the momenta of the fields are transformed according to 




(2.12) 



H(tt, U) = \ (tt, K(U)tt) + S(U) - In det F m (V) 



(2.13) 



where 



K(U) =F.(V) T MV), 



V = F- X {U). 



(2.14) 
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Note that the Jacobian in eq. (2.13) cancels when the momenta are integrated over 
in the functional integral. The algorithm outlined in fig. 1 thus amounts to applying 
the HMC algorithm with a modified Hamilton function of the kind considered long 
ago by Duane et al. [6]. 



3. Transformations generated by flow equations 

Flows in field space build up field transformations from infinitesimal transformations. 
The latter are generally easier to work with than integral transformations, because 
they refer to the current field only. Moreover, the differentiability and invertibility 
of the generated transformations is automatically guaranteed. 

3. 1 Flows in field space 

An infinitesimal field transformation, 

U + eZ(U)U + 0(e 2 ), (3.1) 

is generated by a link field [Z(U)](x, /x) with values in su(3). The continuous com- 
position of such transformations amounts to integrating a flow equation 

U t = Z t (U t )U t (3.2) 

with respect to a ficticious time t. A simple choice for the generator of the flow is 

[Z t (U)r(x,^)=d^W (U), (3.3) 

W (U)= ti{U(x,fi,v)}, (3.4) 

where U (x, fi, v) denotes the plaquette loop in the (n, ^-directions at the point x. 
In the following, this flow will be referred to as the "Wilson flow". 

If Z t (U) is a differentiable function of t and U, the flow equation (3.2) has a unique 
solution U t for any specified initial value Uq = V and all t G (— oo, oo). Moreover, 
the solution is differentiable with respect to t and V (for a proof of these statements, 
see ref. [7] , for example) . It should be emphasized that the existence of the solution 
for all times is non-trivial and can only be guaranteed, without further assumptions, 
because the field manifold is compact. 



6 



3.2 Integrated transformations 

At fixed t, the field U t is a well-defined function of the initial field V. Through the 
integration of the flow equation, one thus obtains a differentiable transformation, 
V — > Ut = J-t(V), of the field space. The transformation is invertible and its inverse 
is differentiable, because the flow equation can be integrated backwards from t to 
0. Moreover, since the Jacobian detJ-'t,*(y) is equal to unity at t = and does not 
pass through zero at any time, the transformation is also orientation-preserving and 
thus fulfils all requirements for an acceptable map of field space. 

There is a useful compact expression for the Jacobian which is obtained starting 
from the equations 



^lndetJT t ,*(F) = Trtft^V^iV)- 1 }, (3.5) 

[^(V)](x, fi; y, ») ab = -2 tr{<9^{[Z t ([/ t )](z, »)U t (x, fi)}U t (x, m)" 1 ^ 

-dlMx^Ptix^y^ZtiUMx,^}. (3.6) 

Noting 

= T.lMV^^iyr^, (3.7) 

a few lines of algebra then lead to the formula 

lndet^(V) = fds 52{dZ tli [Z a (U)] a (x,vi)} u=Ut . (3.8) 



In the case of the Wilson flow, for example, the contribution of the Jacobian to the 
action of the field V, 

\ndetFtAV) = -f [ dsW {U s ), (3.9) 
J o 

is proportional to the integral of the Wilson plaquette action along the flow. 
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4. Trivializing maps 



Somewhat surprisingly, trivializing maps can, to some extent, be constructed expli- 
citly in the pure gauge theory. The construction is explained in this section, assuming 
that the gauge action S(U) is a sum of Wilson loops (plaquettes, rectangles, etc.). 

4-1 Trivializing flows 

If the generator Z t (U) of the flow (3.2) is such that 

fds Y,{d a x AZ s {U)\ a {x^)} u=Ua = tS(U t ) + C t , (4.1) 

where C t may depend on t but not on the fields, the associated integrated transfor- 
mations satisfy 

S(F t {V)) - lndet^ t ,,(F) = (1 - t)S{T t {V)) - C t . (4.2) 

In particular, the transformation at t = 1 is then a trivializing map. 

Equation (4.1) is a rather implicit condition on the generator of the flow. However, 
when differentiated with respect to t, it assumes a more tractable form, 

J2{d a x jZ t (U)] a (x,») -td^S(U)[Z t (U)} a (x,»)} = S(U) + C t , (4.3) 

which involves the generator at time t only. Note that the differential condition (4.3) 
and the flow equation (3.2) imply eq. (4.1), i.e. it suffices to find a generator Z t (U) 
that satisfies eq. (4.3). 

4-2 Existence of trivializing flows 

Equation (4.3) is an inhomogeneous linear partial differential equation for the gener- 
ator Z t (U). Since it is a scalar equation, one expects that there are many solutions. 
In the following, the solution will be obtained in the form 

[Z t (U)] a (x, = J t (U), (4.4) 

where the action S t (U) is to be determined. 
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When inserted in eq. (4.3), the ansatz (4.4) leads to the Laplace equation 



£ t S t = S + C t , (4.5) 



(for simplicity, the argument U is now often omitted). The operator £ t is elliptic 
and symmetric with respect to the scalar product 

(6V) = I B[U]e- ts ^<P(UmU)- (4-7) 

£ t has therefore a complete set of differentiable eigenfunctions and a purely discrete 
spectrum with no accumulation points (see ref. [8], sect. 1.6, for example). Moreover, 
since 

(0, £ t <f>) = ^ ^) > 0, (4.8) 

the function 4>(U) = 1 is the only zero mode of £ t and all other eigenfunctions have 
eigenvalues separated from the origin by a strictly positive spectral gap. 
Now if one chooses C t to be such that 

C t = -(1,S)/(1,1), (4.9) 

the zero-mode component is removed from the right-hand side of eq. (4.5) and 

S t = 2t\S + C t ). (4.10) 

is then a well-defined expression that solves the equation. The differentiability of St 
with respect to t and U essentially follows from the ellipticity of £ t (appendix E). 
A constructive proof of the existence of trivializing flows has thus been given. 

4-3 Expansion in powers oft 

The solution (4.10) is well defined but still quite implicit since it involves the inverse 
of an operator acting on functions of the gauge field. In this and the next subsection, 
the solution is worked out analytically in powers of t. 
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When the series 

oo 

£ t = ^t fc 5 (fc) (4.11) 

fc=0 

is inserted in eq. (4.5), the matching of the terms of a given order in t leads to the 
recursion 

£ S (0) = S + C w , (4.12) 
£ S (k) = J2 d a x ,,Sd a x J^ + C (k \ k = 1, 2, . . . , (4.13) 

for the actions . The Laplacian £ coincides with the colour-electric part of the 
Hamilton operator in lattice gauge theory in 4 + 1 dimensions. In particular, its 
eigenfunctions are products of SU(3) representation functions of the link variables. 
Sums of Wilson loops and products of Wilson loops, for example, are eigenfunctions 
of £o or can easily (algebraically) be decomposed into eigenfunctions. 
The solution of the recursion, 

= £q 1 S, ( 4 - 14 ) 
£(*) = -^Y.d^Sd^J^, k = 1,2,..., (4.15) 

x, p. 

is thus obtained in the form of sums of Wilson loops and products of Wilson loops. 
Note that, as already mentioned in subsect. 4.2, the constant function is the only 
zero mode of £o- Since the smallest non-zero eigenvalue of £o is |, the right-hand 
sides of eqs. (4. 14), (4. 15) are therefore unambiguously determined up to an irrelevant 
additive constant. 

4-4 Calculation of and in the Wilson theory 

For illustration, the first two terms of the series (4.11) are now worked out explicitly 
for the plaquette action [9] 

S W (U) = -\(iW {U) (4.16) 

(where = 6/g$ denotes the inverse gauge coupling). A short calculation, using the 
completeness relation (A. 5), shows that the leading term is 

S<® = -±(3W = ^S w . (4.17) 



10 



(1) 



(2) 



(3) 



(4) 



(5) (6) (7) 

Fig. 2. Classes of loops and pairs of loops contributing to S^ 1 ' in the Wilson theory. 
The loops 5 — 7 reside on a single plaquette of the lattice. All other loops occupy two 
plaquettes which can lie in a plane or be at right angles in three dimensions. 

To this order and up to a rescaling of the time parameter t, the trivializing flow in 
the Wilson theory thus coincides with the Wilson flow. 
The expression to be worked out at the next order is 

S W = - ik^o'E^Wo^Wo- (4.18) 

The Wilson loops and products of Wilson loops that can occur at this point derive 
from the contractions of two plaquette loops with a common link. Altogether there 
are then seven classes Cj, i = 1, . . . , 7, of loops and pairs of loops to consider (see 
fig. 2). By summing the traces of the associated Wilson loops, each class C% defines 
an action 

W t = ^2 tr{U(C)} if * = 1,2, 5, (4.19) 
Wi= Yl tT {U(C)}tr{U(C')} if i = 3,4,6,7, (4.20) 

{C,C'}€C t 

where U (C) denotes the ordered product of the link variables along the loop C. The 
sums in these equations extend over all possible positions of the loops on the lattice. 
Loops with opposite orientation are considered to be different and are both included 
in the sums. 
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Using the identity (A. 5) again, some algebra now yields 

E d l» w ° d ^ w ° = m-m- |w 3 + ±w 4 - 2w 5 + |w 6 - §w 7 (4.21) 

up to an additive constant. Furthermore, 



£ Wl = 8W1, 


(4.22) 


£ W 2 = f w 2 + w 4 , 


(4.23) 


£ W3 = 11W 3 -Wi, 


(4.24) 


£ W 4 = f W 4 + W 2 , 


(4.25) 


£ W 5 = f W 5 + 4W 6 , 


(4.26) 


£ W 6 = f W 6 + 4W 5 , 


(4.27) 


£ W 7 = 12W 7 + constant. 


(4.28) 



In the subspace of these functions, the operator £ can be easily inverted and the 
result 

= Tk? 2 {-&m + mm + j- 3 m - + 

-IW 6 + i>V 7 } (4.29) 

is thus obtained. Note that the smallness of the numerical coefficients in this formula 
is balanced, to some extent, by the number of loops per lattice point in the classes 
Ci (which are equal to 120, 12 and 6 for i = 1, . . . , 4, i = 5, 6 and i = 7 respectively). 

4-5 Miscellaneous remarks 

(a) Higher orders. The actions , k > 2, can be computed algebraically following 
the steps taken in the previous subsection. Since all loops generated by contracting 
a plaquette loop with the loops at order k — 1 must be considered, the work required 
for the calculation is however rapidly growing with k. 

(b) Locality and convergence. The series (4.11) is an expansion in local terms whose 
footprint on the lattice increases proportionally to the order k. At the values of t, 
where the expansion converges, the action S t is then guaranteed to be local as well. 
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The norm estimates in appendix E imply a lower bound on the convergence radius 
of the series, but this bound is rather poor and vanishes in the infinite-volume limit. 
It seems nevertheless plausible that the series has a non-zero convergence radius in 
this limit, because the inverse of the operator £< in eq. (4.10) is likely to remain 
bounded in a complex neighbourhood of t = 0. An analysis that takes the locality 
properties of £ t into account will however be required to show this. 

(c) Truncation of the expansion. If all terms in the series (4.11) of order k > n are 
dropped, one obtains an approximately trivializing flow that satisfies eq. (4.2) up to 
an additive correction of order t n+1 . An at least partial cancellation of the action is 
achieved in this case as long as t is sufficiently small for the correction to be strongly 
suppressed. 

(d) Smoothing property. The Wilson flow satisfies 

^S w (U t ) = -± ^{^S w (U)^S w (U)} u=Ut < (4.30) 

and therefore lowers the Wilson action as t increases. To leading order in t, the triv- 
ializing flow constructed in this section for the Wilson theory thus has a smoothing 
effect on the gauge field. On the other hand, if the flow is followed in the reverse 
direction, the gauge field tends to become rougher. 

(e) Topological charge sectors. In lattice QCD, the topological (instanton) sectors are 
not a property of the field manifold alone, but are expected to emerge dynamically 
when the continuum limit is approached. The fact that trivializing maps completely 
"straighten out" the sectors is therefore not in conflict with the topological properties 
of the field space. 

(f) Renormalization group. By composing the trivializing map U = ^i(V) in the 
Wilson theory with its inverse at another value of the gauge coupling, one obtains a 
group of transformations whose only effect on the action is a shift of the coupling. 
The locality properties of these transformations are not transparent, however, and 
could be quite different from the ones of a Wilsonian "block spin" transformation. 



5. Numerical integration of the Wilson flow 

The discussion in sect. 1 now suggests to combine the HMC algorithm with the field 
transformations generated by the trivializing flow constructed in the previous sec- 



13 



tion. In particular, if the Wilson gauge action is used, the transformations generated 
by the Wilson flow may lead to an algorithm with improved sampling efficiency 



5.1 Forward integration 

There is a wide range of numerical integration methods that can in principle be used 
to integrate the Wilson flow (see ref. [10], for example). The Euler scheme discussed 
in the following performs the integration in time steps of size e and updates the link 
variables one after another according to 

U(x,n) -» U'{x,n) = e e ^ z ™ x ^U(x,n), (5.1) 

where 

[Z(U)](x,ri=T a dl tl Wo(U). (5.2) 

Starting from the gauge field U t at time t, the field at time t + e is thus obtained by 
running through all links (x, /j.) on the lattice and updating the link variable residing 
there. Note that the ordering of the links matters, since the old value of U(x,fi) is 
replaced by the new one before going to the next link. 
The generator of the flow is explicitly given by 

[Z(U)](x, fi) = -J2 V{U{x, n)U(x + A, v)U{x + z>, fi^Uix, v)~ x 

+ U(x,fj,)U(x + /t - i>, v^Uix - u,fi)~ 1 U(x -v,v)}, (5.3) 

where (i denotes the unit vector in direction \i and 

V{M) = \ (M - M 1 ") - £tr(M - M r ) (5.4) 

projects any 3x3 matrix M to su(3). The Euler integration of the Wilson flow thus 
amounts to applying a number of "stout smearing" steps [11], except that the link 
variables are here updated one by one rather than all at once. 

5.2 Backward integration 

The application of n Euler sweeps maps the initial field V = Uq to the field U = U ne 
at time t = ne. If t is held fixed and n is taken to infinity, this map converges to the 
transformation obtained by integrating the Wilson flow exactly. However, the HMC 
algorithm may potentially be combined with the map defined by the Euler integrator 
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at fixed n and e. For this proposition to be a viable option, the transformation must 
be invertible, i.e. one must be able to trace back the Euler integration by inverting 
the link update steps one by one in the reverse order. 

The question is thus whether eq. (5.1) has a unique solution U(x, fx) given U'(x, fx) 
(and keeping all other link variables fixed). As explained in appendix D, the answer 
is affirmative, for arbitrary values of the field variables, if 

|e| < I (5.5) 

Moreover, in this range of e, the solution U{x,n) = e~ eX *U'(x, fi) can be obtained 
through the fixed-point iteration 

X Q = 0, (5.6) 

X n+1 = {[Z(U)](x,»)} UM=e _ eXnU , M , n = 0,1,2,..., (5.7) 

which converges to X* at an exponential rate. 

In appendix D it is also shown that the Jacobian of the transformation (5.1) is 
strictly positive in the range (5.5). The field transformations obtained through the 
Euler integration of the Wilson flow are therefore orientation-preserving diffeomor- 
phisms of the field manifold, as required for acceptable maps of field space. 

5.3 Jacobian matrix of the Euler integrator 

The Euler step (5.1) amounts to applying the field transformation 

f e^ z( - u ^ x ^U(x, n) if (y, v) = (x, ft), 

[ U (y, v) otherwise, 

to the current gauge field. An Euler sweep is then the composition product of these 
transformations over all links (x,/i). It is straightforward to show that the Jacobian 
matrix of a composed transformation is the product of the Jacobian matrices of the 
factors. The Jacobian matrix of the Euler integrator is therefore an ordered product 
of the Jacobian matrices 

[£ w (U)}(y^;z,pr = -2tT{dl p [£ x ^{U)]{y,u)[£ x ^{U)]{y,u)- l T a } (5.9) 

of the one-link transformations (5.8). 

The matrix (5.9) can be expressed through the derivative 

[Z*(U)](y,v;z,p) bc = d c z jZ(U)] b (y,u) (5.10) 
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of the generator of the Wilson flow. Explicitly, one finds that 



[£ x ,^(U)}(y^;z,p) ac = 5 ac 5 yz 5 up + 5 xy 5^[{e AdX - l) a ° 5 yz 5 vp 

+ eJ{-XY b [Z*{U)]{y,v-z,p) bc \ (5.11) 

> X=e[Z(U)\(x,fi) 

where use was made of the SU(3) formulae listed in appendix A and B. 

5.4 Jacobian of the integrated transformations 
The Euler integrator generates a sequence of fields 

V = U -> U e -> U 2e -> . . . - U ne = U (5.12) 

by sweeping through the lattice n times and updating the link variables one by one 
in a specified order. In the following, the intermediate field configurations obtained 
starting from Uk e and updating the link variables on all links that come before (x, fi) 
will be denoted by U ke ^ X fl j. In particular, U ke ^ X: ^ = Uke if (x,n) is the first link 
and 



Uke,[y,v] = £xA U ke,[x,iJ,]) (5-13) 

if (y, v) follows (x, jj) in the chosen link order. 

Since the transformation V — > U = J- n eiy) is a composition product of one-link 
update steps, its Jacobian 

n-l 

detTneAV) = II H det£ ^*(Uke,M) (5-14) 

k=0 x,fj, 

factorizes into the product of the Jacobians of the steps. The latter coincide with 
the determinants of certain real 8x8 matrices given explicitly in appendix C. 



6. Proposed simulation algorithm for lattice QCD 

With respect to the QCD simulation algorithms used to date, the combination of the 
transformations obtained through the Euler integration of the Wilson flow and the 
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HMC algorithm is expected to sample the topological sectors more quickly and to be 
generally more efficient. The use of the Wilson flow is suggested if the gauge action 
coincides with the Wilson plaquette action. For other actions, the appropriate flow 
can be easily constructed following the lines of sect. 4. 

6. 1 Choice of the parameters 

The parameters of the Euler integrator are the integration step size e and the number 
n of Euler sweeps that are applied. One also needs to choose a definite ordering of 
the links of the lattice. 

The step size e should be positive and not larger than, say, 1/16 so that the inver- 
tibility of the Euler integrator is guaranteed within a safe margin. Some tuning of 
the integration time ne will certainly be required in order to maximize the efficiency 
of the algorithm. Note that the unit of time differs from the one used in subsect. 4.4, 
i.e. setting t = 1 there corresponds to an integration time ne = (3/2>2. 

The ordering of the links can be chosen arbitrarily. One may, for example, first 
visit the links (x, 0) on all even points x, then the links (x, 0) on all odd points, then 
the links (x, 1) on all even points, and so on. This ordering is well suited for parallel 
processing, since the link variables in a given direction on the even (odd) sites are 
decoupled from one another and can therefore be updated in parallel. 

6.2 Force calculation 

At the beginning of an update cycle, the current gauge field U is transformed to the 
field V = J-~^{U) by applying n backward Euler sweeps to U. The force that drives 
the molecular-dynamics evolution of V is then given by 



where the action S(U) now includes the usual sea-quark pseudo-fermion actions (for 
simplicity, the dependence on the pseudo-fermion fields is suppressed). 

Each time the force is to be calculated, the current field V must be transformed to 
U again by applying n forward Euler sweeps (see fig. 3) . The fields Uq,U £ ,..., U ne 
generated in this process should be stored in memory so that they will be available 
when the force is propagated from U to V. Note that the intermediate fields 



F(z,p) c = dl p {S(T ne (V)) - lndetJ^F)} , 



(6.1) 




(6.2) 



are then also available. 
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Fig. 3. The proposed algorithm evolves the field V using the standard HMC algo- 
rithm (thick line) . The force that drives the molecular-dynamics evolution is obtained 
by forward integration of the Wilson flow and subsequent backward propagation of 
the derivatives of the action and the Jacobian of the flow (n = 3 in this figure) . 

The factorization (5.14) of the Jacobian implies 

n-l 

F(z, P y = Bl p S(U ne ) - ^^^lndet^^,^]). (6.3) 

k=0 x.fi 

Moreover, 

d c z>p S(U ne ) = dl v S{U)\ u=Un T nej *{y,v;z,p) ac , (6.4) 

and there is a similar formula for the other terms in eq. (6.3) involving the Jacobian 
matrices of the transformations V — > Uke,[x,n]' Ah 1 these matrices, as well as the one 
in eq. (6.4), are products of the Jacobian matrices of the one-link transformations 
(5.8). The force can therefore be computed recursively as follows: 

1. Set U = U ne and F(z,p) c = d c Zjf) S(U). 

2. For k from n — 1 to 0, run through all links (x, jx) in reverse order, set U = Uk e ,[ x ,^] 
and update the force according to 

Note that the field U backtracks the forward integration of the Wilson flow in the 
course of the recursion. Since the Jacobian matrix £ x ,u,*(U) differs from unity only 
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on the links sharing a plaquette with the total computational effort required 

for the propagation (6.5) of the force is expected to be similar to the one required 
for n applications of a nearest-neighbour gauge-covariant difference operator to the 
force field. 

6.3 Domain-decomposed algorithm 

Field transformations can also be combined with the DD-HMC algorithm [12]. Only 
the so-called active link variables are transformed in this case, but the transforma- 
tions may depend on the inactive field components. 

In the pure gauge theory, it is then again possible to construct trivializing flows. 
They operate on the active link variables and contract the gauge action in each do- 
main to a constant (i.e. an expression depending on the inactive link variables only). 
These flows are not the same as the ones constructed in sect. 4, but can be obtained 
by solving the differential equations derived there. In particular, the trivializing flow 
for the Wilson action is, to leading order, a slightly modified Wilson flow, where the 
plaquettes containing v active links are given the weight 4/za 



7. Concluding remarks 

Trivializing maps and flows in field space have here been discussed with a particular 
application in mind. The underlying concepts are fairly general, however, and may 
have other uses in rigorous constructive work, numerical perturbation theory or in 
connection with renormalizable smoothing techniques, for example. 

Whether the proposed combination of the Wilson flow and the HMC algorithm 
does indeed sample the topological sectors in lattice QCD more efficiently than the 
simulation algorithms used so far remains to be determined. As the lattice spacing 
is taken to smaller and smaller values, the quark fields may eventually have to be 
included in the flow and perhaps also the next-to-leading order correction discussed 
in sect. 4. Trivializing maps in the presence of matter fields is, in any case, a subject 
that deserves to be studied in its own right. 

In the course of this work, I profited from many discussions with Filippo Palombi 
and Stefan Schaefer of various questions related to the slow topology-switching in 
current lattice QCD simulations. I also wish to thank Stefan Schaefer, Rainer Som- 
mer and Francesco Virotta for sharing some of their simulation results before pub- 
lication. 
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Appendix A. SU(3) notation 



A.l Group generators 

The Lie algebra su(3) of SU(3) may be identified with the space of all anti-hermitian 
traceless 3x3 matrices. With respect to a basis T a , a = 1, . . . , 8, of such matrices, 
the elements X £ su(3) are given by 

X = X a T a , (A.l) 

where (X 1 , . . . ,X 8 ) £ M 8 (repeated group indices are automatically summed over). 
The generators T a are assumed to satisfy the normalization condition 

tr{T a T b } = -\5 ab . (A.2) 
The structure of the Lie algebra is then encoded in the commutators 

[T a ,T b ] = f abc T c , (A.3) 
while the completeness of the generators implies 

{T a , T b } = -±5 ab + id abc T c , (A.4) 

T a(S T yS = -§{<5a5<*/37 ~ k^Sjs} ■ (A.5) 

It follows from these equations that the structure constants f abc and the tensor d abc 
are both real. Moreover, f abc is totally anti-symmetric in the indices and d abc totally 
symmetric and traceless. 

A.2 Adjoint representation 

The representation space of the adjoint representation of su(3) is the Lie algebra 
itself, i.e. the elements X of su(3) are represented by linear transformations 

AdX : 5u(3) ^su(3), (A.6) 

AdX-Y = [X,Y] for all Y € su(3). (A. 7) 

The action of Ad X on the group generators is given by 

AdX -T b = T a (Ad X) ab , (A.8) 
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where 

(AdX) ab = -f abc X c (A.9) 
is a real antisymmetric 8x8 matrix. 
A. 3 Matrix norms 

The natural scalar product in su(3) is 

(A, Y) = X a Y a = -2 tr{XY}. (A.10) 

In particular, ||A| = (X, X) 1 / 2 is a possible definition of the norm of X € su(3). 
Another useful matrix norm derives from the square norm 

IM| 2 = {bi| 2 + M 2 + N 2 } 1/2 (A.ii) 

of complex colour vectors v. If A is any complex 3x3 matrix, one defines 

P|| 2 = max ||Au|| 2 . (A.12) 

1 1 U 1 1 2 = 1 

This norm satisfies 

\\A + B\\ 2 < \\A\\ 2 + \\B\\ 2 , \\AB\\ 2 < \\A\\ 2 \\B\\ 2 , (A.13) 



for all matrices A, B. Moreover, if A is hermitian or antihermitian, ||^4||2 is equal to 
the maximum of the absolute values of its eigenvalues. 



Appendix B. Properties of the SU(3) exponential function 

B.l Lipschitz bound 

For any 1,7 6 su(3), the relation 

11^-^2 = ||1 -e-*e Y || 2 (B.l) 

follows from the fact that e x is unitary. Using the identity 

1 - e~ x e Y = f ds e~ sX (X - Y)e sY (B.2) 
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and the subadditivity (A. 13) of the norm, the Lipschitz bound 

||e x -e y || 2 < ||A-Y|| 2 (B.3) 
is then obtained. 

B. 2 Differential of the exponential map 

Let X be an element of the Lie algebra su(3). A linear mapping J(X) : su(3) i— > su(3) 
is then defined by 

for all Y € su(3). (B.4) 



J(X)-Y = e~ x ±e x+tY 
at 



t=o 



J(X) is referred to as the differential of the exponential map. Scaling the exponents 
by a parameter s, as above, one obtains the representation 



J(X)-Y = [ dse~ sX Ye sX 
Jo 



(B.5) 



and thus the expansion 



00 (—-\\ k 

J(X) = l + ^f T ^(^X) k (B.6) 

which is absolutely convergent for any X £ su(3). 

The action of J(X) on the group generators T a is given by 

J(X) -T b = T a J(X) ab , (B.7) 
where J(X) ab is a real 8x8 matrix. Note that 

J(X) T = J(-X) = e AdX J(X), J(X)AdX = l-e~ AdX . (B.8) 
Moreover, eq. (B.5) implies 

\\J(X) -Y\\ 2 < \\Y\\ 2 (B.9) 
for all X,Y G su(3). 
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Appendix C. Jacobian of the Euler step 

The Jacobian matrix (5.11) of the Euler step is equal to unity except for some non- 
zero elements along the row (y, v) = (x,n). Its determinant therefore coincides with 
the determinant of the (x, fi; x, /i)-element 

A ac = \(e Adx ) ac + eJ(-X) ab [Z*(U)](x, t i;x,Li) bc } (C.l) 
of the matrix. 

The derivative [Z*{U)](x, /i; x, n) bc can be worked out explicitly in terms of the 
plaquette sum 

M = ^2{U(x, n)U(x + A, v)U(x + v, fi)~ 1 U(x, + 

U(x, fi)U{x + /t - z>, v)~ l U{x - z>, n)~*U{x -i>,v)}. (C.2) 

A few lines of algebra then lead to the expression 

A ac = B ac + i eC ab {id bcd tr{T d (M + M f )} - \5 bc tx{M + M f }} (C.3) 
where 



B = l(e AdX + l), 

C = J(-X), X = -eV{M}. 
In particular, 

det A = 1 - |etr{M + M 1 "} + 0(e 2 ), 
as expected from eq. (3.9). 



(C.4) 
(C.5) 

(C.6) 
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Appendix D. Inversion of the Euler step 

D.l Basic norm bounds 

For any complex 3x3 matrix M, the inequality 

\\V{M}\\ 2 < l\\M\\ 2 (D.l) 
can be established in a few lines. One first observes that 

i(M-M ! ) =ADA~\ i G SU(3), (D.2) 

where D is a diagonal matrix with diagonal elements Ai, A2, A3. Setting 

3 

A=§X>, (D.3) 

k=l 

the estimates 

||P{M}|| 2 =max|A fc -A| 

k 

< |max|A fc | = -Aft|| 2 < |||M|| 2 (D.4) 

then show that the inequality (D.l) holds for all matrices M. 

An immediate consequence of the bound (D.l) and the Lipschitz bound (B.3) is 
that 

\\V{A(e x -e Y )B}\\ 2 <l\\X-Y\\ 2 (D.5) 
for all A,B e SU(3) and 1,7 6 su(3). 
D.2 Solution of eq. (5.1) 

For a given link (x,fi) and any fixed U'(x,fi) <G SU(3), the function 

f(X) = {[Z(U)](x,»)} UM=e _ eXu , M (D.6) 

maps X G su(3) back to su(3). Recalling eq. (5.3), the inequality (D.5) immediately 
implies that 

||/(X)-/(y)|| 2 <A;||X-y|| 2 , fc = 8|e|. (D.7) 

24 



The function f(X) is therefore a strict contraction if the integration step size e is in 
the range (5.5) (which is assumed to be the case from now on). 

It is not difficult to prove that strict contractions in a complete metric space have 
a unique fixed point (see ref. [13], Theorem V.18, for example). In the present case, 
the fixed point X. ¥ can be computed by noting that the sequence X = 0, Xx, X 2 , . . . 
generated through the recursion X n+1 = f(X n ) satisfies 

\\X n - X*\\ 2 < k\\X n -! - X*\\ 2 CD- 8) 

and therefore rapidly converges to X*. The matrix 

U(x,n) = e- eX *U'(x,n) (D.9) 

then provides a solution of eq. (5.1). Moreover, there is no other solution, because 
the fixed point X* of f(X) is unique. 

D.3 Positivity of the Jacobian 

In order to prove that the determinant of the Jacobian matrix (C.l) is positive at all 
e in the range (5.5), it suffices to show that the matrix has no zero mode, i.e. that 

A ac Y c = (D.10) 
implies Y = 0. To this end, eq. (D.10) is first written in the form 

Y a = -eJ(X) ab W b , (D.ll) 
where 

X = e[Z(U)](x,fi), W b = [Z*(U)](x, f i- 1 x, f i) bc Y c . (D.12) 
Recalling eq. (5.10), the formula 

W = Si { W U N x >ri\u(^)^u M - [Z(U)Kx,f)} (D.13) 

may then be derived from which one infers that 

||W|| 2 <8||Y|| 2 . (D.14) 

The inequality (D.5) has here been used again and also the fact that the norm is 
a continuous map from su(3) to R. The combination of eq. (D.ll) and the bounds 
(B.9) and (D.14) now implies ||Y|| 2 < fc||Y|| 2 and thus Y = 0. 
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Appendix E. Differentiability of S t 



The action S t solves an elliptic partial differential equation with smooth coefficients 
and is therefore guaranteed (by elliptic regularity) to be a differentiable function of 
the gauge field U. In this appendix, the simultaneous differentiability in the time t 
is established by expanding S t in powers of t — t around any fixed time t . Using 
Sobolev norms, the expansion can be shown to converge if t is sufficiently close to to. 
The pointwise uniform convergence of the series and its derivatives, and therefore 
the differentiability of St, then follows from Sobolev's lemma. 

E.l Sobolev spaces on the field manifold 

The definition of the Sobolev spaces on a compact manifold is quite involved and will 
not be reviewed here. An introduction to the subject and the proof of all statements 
made in this subsection is given in the first three sections of ref. [8], for example. 

Let C°° be the space of differentiable functions on the field manifold and the 
associated Sobolev space of order k £ Z. The latter is the completion of C°° with 
respect to a certain norm || • ||fe. A characteristic feature of these norms is that the 
bounds 



Wh < c k U\\ k+p (E.l) 

hold for all (f> £ C°° and any differential operator D of order p with smooth coeffi- 
cients, where the constants c k depend on D but not on <p. Such differential operators 
thus extend to bounded linear operators from Hk +P to H^. Moreover, \\4>\\k < \\(t>\\i 
if k < I and therefore Hi C Hh- 

A fairly concrete description of the Sobolev space Hk can be given when k is even 
and non-negative. For any tel and j = 0, 1, 2, . . ., another norm \\4>\\t,j of G C°° 
may be defined through 

U\\t,j = W\ + (E.2) 

where || • || is the norm associated to the scalar product (4.7). The fact that is a 
second-order elliptic differential operator then implies 

IH| 2j <a tj IHb', UWtj^hjUhj (E.3) 

for some constants a t j, b tt j. In particular, H 2 j is the completion of C°° with respect 
to the norm || • \\ t j. 
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E.2 Properties of the inverse of £ t 

Let tyt be the orthogonal projector to the zero mode of £ t in the Hilbert space with 
scalar product (4.7). Its action on any function cf) <G C°° is given by- 



Note that <p t projects to a constant, but the constant depends on the definition 
of the scalar product (• , •) and therefore on t. The action of the inverse &t of £t on 
(f> is then determined by the equations 



As already mentioned in subsect. 4.2, the ellipticity of £t and the absence of further 
zero modes imply that these equations have, for any fixed i, a unique differentiable 
solution Evidently, S t = ^tS. 

Since £ t has a spectral gap in the subspace orthogonal to the zero mode, & t is a 
bounded operator with respect to the norm || • ||. As a consequence, there exists a 
constant gt such that 



for all j = 1,2, . . . and all <j) £ C°°. Moreover, recalling the inequalities (E.3), one 
concludes that 



= </>)/(!, 1). 



(E.4) 



£ t «S t = (1 - ^ t )<f>, 



y t ®t<t> = o. 



(E.5) 



®t<t>\\t,j < 9tU\\t,j-i 



(E.6) 



®t<t>hj < gtjUhj-2 



(E.7) 



for some other constants g t j . 



E. 3 Expansion of S t 

For any fixed time to, the operator £ t may be decomposed according to 





(E.8) 



Each term in the power series 



oo 




(E.9) 



n=0 
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is then a well-defined differentiable function of both t and U. Moreover, since 

II W>lk 3 < gtojlWhj-2 < djgtojthj-i < djgt^Hhj, (E.10) 

all j = 1, 2, . . . and some constants dj, it is clear that the series and all its derivatives 
with respect to t converge uniformly in the norm || • if i is sufficiently close to 
to- Sobolev's lemma (statement (c) in lemma 1.3.5 of ref. [8]) then implies that the 
series converges pointwise and uniformly, together with all its derivatives in t and 
its derivatives in U up to some order proportional to j. 

If t is in a neighbourhood of to, where the convergence of the series and its deriva- 
tives up to order I > 2 is guaranteed, the action of the operator £ t and the summation 
may be interchanged and one finds that 

oo 

£tVt = S - ~ tf <P* (£'<S to ) n S = (1 - «p t ) S, (E.ll) 

n=0 

the second equality being implied by the identities y$ t £, t ip t = and tyttyt = tyt - 
As a consequence, 

{\-%)^ t = & t S = S u (E.12) 

which shows that St is, in the specified neighbourhood of to, I times continuously 
differentiable with respect to t and U. Since to and I can be chosen arbitrarily, the 
simultaneous differentiability of St is thus guaranteed at all times t and to all orders. 
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